A8Y 

■NICAL  REPORT  SECTfO* 

AVAL  POSTGRADUATE  JCWnn, 
-O^R^CALXTO^^0' 


NPS-69BC76101 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


FINITE  DIFFERENCE  TREATMENT  OF  TRANSIENT 
TEMPERATURE  DISTRIBUTION  IN  AN  INSULATED  PIPE 


by 


John  E.  Brock 


October  1976 


FEDDOCS 

D  208.14/2:NPS-69BC76101 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  Isham  Linder  T  D  n 

Superintendent  J-   R'  Borsting 

Provost 


FINITE  DIFFERENCE  TREATMENT  OF  TRANSIENT 
TEMPERATURE  DISTRIBUTION  IN  AN  INSULATED  PIPE 

This  report  discusses  a  Crank-Nicolson  type  of  finite 
difference  formulation  for  the  problem  of  transient 
temperature  distribution  in  a  uniform  cylindrical  pipe 
with  external  insulation  and  containing  a  fluid  for 
which  temperature,  pressure,  and  flow  rate  are  given 
as  functions  of  time.  All  material  properties  are 
assumed  to  depend  on  temperature.  Discussion  is  also 
given  of  a  digital  computer  program  which  implements 
this  formulation.  This  implementation  incorporates 
properties  of  superheated  steam,  the  internal  fluid, 
and  of  2-1/4  Cr  1  Mo  (P22),  the  pipe  material.   Stress 
calculations  are  also  treated. 


NPS-69BC76101 
October  1976 


UNCLASSIFIED 


SECURITY   CLASSIFICATION   OF   THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1.     REPORT   NUMBER 

NPS-69BC76101 

2.  GOVT   ACCESSION  NO. 

3.     RECIPIENT'S  CATALOG   NUMBER 

4.     TITLE  (and  Subtitle) 

FINITE  DIFFERENCE  TREATMENT  OF  TRANSIENT 
TEMPERATURE  DISTRIBUTION  IN  AN  INSULATED  PIPE 

5.     TYPE  OF   REPORT  ft   PERIOD  COVERED 

6.     PERFORMING  ORG.   REPORT  NUMBER 

7.     AUTHORS 

John  E.   Brock 

8.     CONTRACT  OR  GRANT  NUMBERfs) 

9.     PERFORMING  ORGANIZATION   NAME  AND  ADDRESS 

Professor  John  E.   Brock  (Code  69Bc) 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 

10.     PROGRAM   ELEMENT,  PROJECT,   TASK 
AREA  ft  WORK  UNIT  NUMBERS 

11.     CONTROLLING  OFFICE  NAME   AND  ADDRESS 

12.     REPORT   DATE 

October  19 7 6 

13.     NUMBER  OF  PAGES 

U-     MONITORING  AGENCY  NAME  ft    ADDRESSf//  different  from  Controlling  Olfice) 

15.     SECURITY  CLASS,  (of  this  report) 

Unclassified 

15a.     DECLASSIFI  CATION/ DOWN  GRADING 
SCHEDULE 

16.     DISTRIBUTION   ST  ATEMEN  T  (of  this  Report) 

Approved  for  public  release;  distribution  unlimited. 

17.     DISTRIBUTION  STATEMENT  (of  the  abstract  entered  In  Block  20,  if  different  from  Report) 

18.     SUPPLEMENTARY   NOTES 

19.     KEY  WORDS  (Continue  on  reverse  aide  if  necessary  and  Identity  by  block  number) 

Transient  temperatures  in  pipe 
Transient  radial  temperature 
Nonlinear  temperature  transients 

20.     ABSTRACT  (Continue  on  reverse  side  If  necessary  and  Identify  by  block  number) 

This  report  discusses  a  Crank-Nicolson  type  of  finite  difference  formulation 
for  the  problem  of  transient  temperature  distribution  in  a  uniform  cylindrical 
pipe  with  external  insulation  and  containing  a  fluid  for  which  temperature, 
pressure,  and  flow  rate  are  given  as  functions  of  time.     All  material  proper- 
ties are  assumed  to  depend  on  temperature.     Discussion  is  also  given  of  a  digi- 
tal computer  program  which  implements  this  formulation.     This  implementation 
incorporates  properties  of  superheated  steam,  the  internal  fluid,  and  of  2-1/4 
Cr  1  Mo  (P22),  the  pipe  material.      Stress  ralculati nns  app  r!p.c>  f-pp.atpd. 

DD     |   jan"73     1473  EDITION  OF    1   NOV  65  IS  OBSOLETE 

S/N   0102-014-  6601   | 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Bntarad) 


FINITE  DIFFERENCE  TREATMENT  OF  TRANSIENT  TEMPERATURE 
DISTRIBUTION  IN  AN  INSULATED  PIPE 

by  J.  E.  Brock 

Abstract:     This  report  discusses  a  Crank-Nicolson  type  of  finite 
difference  formulation  for  the  problem  of  transient  temperature 
distribution  in  a  uniform  cylindrical  pipe  with  external  insulation 
and  containing  a  fluid  for  which  temperature,  pressure,  and  flow 
rate  are  given  as  functions  of  time.     All  material  nronerties  are 
assumed  to  depend  on  temperature.     Discussion  is  also  given  of  a 
digital  computer  program  which  implements  this   formulation.     This 
implementation  incorporates  properties  of  superheated  steam,  the 
internal  fluid,  and  of  2-1/4  Cr  1  Mo  (P22),  the  pipe  material. 
Stress  calculations  are  also  treated. 
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1.  Objective 

The  purpose  of  this  monograph  is  to  present  and  to  document  a 
procedure  for  determining  transient  temperature  distribution  and  the 
corresponding  significant  thermal  stress  components  in  metallic  pipes 
which  are  insulated  externally  and  v/hich  contain  steam  the  tempera- 
ture, pressure,  and  flow  history  of  which  is  specified.     Although  a 
classical  solution  is  known  (13)*  for  an  idealized  version  of  this 
problem,  it  is  not  capable  of  dealing  with  the  variable  heat  trans- 
fer coefficient  which  governs  the  exchange  of  heat  between  the  pipe 
wall  and  the  flowing  steam  it  contains;  this  coefficient  is  a  comp- 
licated function  of  the  properties  of  steam.     Furthermore,  the  an- 
alytic solution  presumes  constant  metal  thermal  conductivity  and 
thermal  diffusivity  while  actual  piping  materials  exhibit  very  sig- 
nificant variations  of  conductivity  and  diffusivity  over  ranges  of 
temperature  usually  encountered  in  engineering  practice.     Accordingly 
a  numerical  approach,  capable  of  dealing  with  these  variations,  has 
been  chosen  for  the  procedure. 

2.  Derivation  of  the  governing  nonlinear  partial  differential  equation 

A  convenient  point  of  departure  for  the  derivation  of  the  govern- 
ing nonlinear  partial  differential  equation  is  the  relation     (  5  ) 

,  9T       9   ,uy,  91\        9    ,vX.  9TN        9    ,Xu,  9T\  ,,x 

X^pc9t  =  I^rt?5  f  9^9^  +  9^9^  (1) 

Here  and  later  T  denotes  temperature,  t  denotes  time,  c  denotes  spec- 


*  References  are  listed  on  pages  28  and  29. 
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ific     heat  at  constant  pressure,  and  p  denotes  mass  per  unit  volume. 
Here  £,  n,  and  c  denote  generalized  space  coordinates  and  A,  w,  .and 
v  are  coefficients  such  that  the  distance  ds  between  points   (£,n,C) 
and  (F,+d£,n+dn,C+dO     is  given  by 

(ds)2=  A2(d02  +  y2(dn)2  +  v2(drj2  (2) 

Specializing  to  the  cylindrical  coordinate  system,  we  take  C  =  r, 
n  =  0,   z,  =  z.     Furthermore ,  we  assume  axial  symmetry  so  that  ~  =  0 

a  " 

and  we  presume  no  variation  with  respect  to  axial  position  so  that 

ST 

jZ  =  0.     We  thus  easily  arrive  at  the  evaluations  X  =  v  =  1,  p  =  r, 

and  at  the  equation 

1  9T  =  97T         1  9T       1  9k  9T  ,. 

a  3t       9T2"        r  3r       k  dr  Br  {i) 

where  a  =  k/pc  is  the  diffusivity.     The  variation  of  k  with  resoect 
to  r  is  due  solely  to  the  fact  that  conductivity  is  a  function  of 
temperature  which  itself  depends  on  radius.     Thus  we  finally  arrive 

^  I2£.^  +   1  il+   (k\(9T)2  {Li) 

a  3t       3r2       r  3r  Mk  M9r;  ^  J 

where  k'   =  dk/dT  for  the  pipe  material.     If  the  variation  of  k  with 

respect  to  T  were  simple,  there  would  be  some  point  in  dealing  with 

the  Kirchhoff  transformation  (6).     However,  it  is  not,  and  we  will 

deal  directly  with  the  nonlinear  partial  differential  equation  'I. 

3.  Finite  difference  formulation 

We  employ  a  Crank-Nicolson  finite  difference  formulation,   as 

follows.     Equation  4  is  differenced  in  time  according  to 

[(UT)+  +   (UT)]/2  =   (T+  -  T)/[6(a+  +  a)/2]  (5) 

where  6  represents  the  time  increment,  superscript  (  )  Indicates 

evaluation  at  the  advanced  time  t  +  6,  and  U  is  the  space  operation 
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The  latter  is  next  represented  approximately  by  a  finite  difference 
operation.     Vfe  use  the  notations 


A  =  radius  increment,  A  =  — =- 

8rz 


,  p-  k'A 


ri 


n  _  9T 


(Note  the  new  use  of  symbol  p.)  Then,  generally, 

T.,,  =  T.  +  DA  +  AA2/2;  T.  .  =  T,  -BA  +  AA2/2  (7) 

l+l    l  '   l-l    l 

where  subscript  (.)  indicates  evaluation  at  r  =  r.  =  n.  .  subscript  (.  , ) 

l  i    i  *  l-l 

indicates  evaluation  at  r  =  r.  ,  =  (n.-l)  ,  etc.  (In  general  n.  will  not 
be  an  integer. )  We  easily  find 


A  =  (T±+1  +  Ti_1  -2Ti)/A2;   B  =  (T±+1  -  T._1)/2A       (8) 


so  that 


HT  =  (T^^.^T.  ♦  (^^..^[1/2^  +  PiWyrt.j)/*]}/*1    (9) 
Setting  3.  =  i4A2/6(a.+a. ) ,  we  arrange  the  general  equation  in  the  form 

l 

=  (3.-2)T.  +  T,...  +  T.  ,  +  (T.^-T,  _)[JL  +  p.(T.A,-T,  ,)/'!]  (10) 
l    i    l+l    l-l     l+l  i-l/L2n.   Ki  l+l  i-1"  J 

On  the  left,  the  bracketed  coefficients  of  the  advanced  temperatures 

+       + 

T.  1  and  T. .,  themselves  contain  these  (unknown)  advanced  temperatures. 

Accordingly,  it  is  contemplated  that  iteration  must  be  used  in  which  earlier 
evaluations  of  T.  ,  and  T.  ,  are  used  in  evaluating  the  bracketed  coeffi- 
cients . 

Equation  10  is  valid  for  interior  nodal  points,  i  =  2,3,...,N,  where 
N  is  the  number  of  subdivisions  of  the  pipe  wall.  At  the  inside  surface 
where  the  radius  is  r,  =  n,A,  there  is  a  convective  boundary  condition 


ST 
k^  =  h(T1-4)  (ii) 

where  <$>  is  the  fluid  tenperature  and  h  is  the  surface  heat  transfer 
coefficient.     We  write  this  as 

=  B  =   (a/A)(T  -d>)  (12) 


9T 
9r 


where  a  =  hA/k. ,  and  we  also  use  the  first  of  equations  7,  taking 
i  =  1.  Thus  we  get 

UT  =  {2(Ta-Tj)  -  a(Tr<D)[2  -  i-  -  D^IT,^)]}^        (13) 
Thus,  the  appropriate  equation  for  the  inner  surface  is 
(31  +  2  +  o+[2  -  jL  _  ota+(Tt-2(J)+)]}Tt  -  2tJ 

=  (31-2)T1  +  2T2  -  o(Tx-(j))[2  -  £-  -  o^T  -<!>)] 

+  o+4>+(2   -i-+  pVd»+)  (111) 

i 

where  3X  =  24A2/5(a1+a  ).     Again,  note  the  nonlinearity  represented 
by  the  appearance  of  the  advanced  temperatures  in  the  bracketed  coef- 
ficients and  on  the  right  side.     The  advanced  tenperature  also  appears 
in  B1  and  in  a     since  otj     and  h  depend  on  the  advanced  inside  surface 
metal  temperature. 

Finally,  the  last  equation  pertains  to  the  exterior  surface  node 

9T 
(numbered  M  =  N+l)  where  there  is  perfect  insulation  so  that  ■*—  =  0. 

Employing  also  the  second  o£  equations  7  with  i  =  M,  we  get 

HT  =  2(TN-TM)/A2  (15) 

and  the  last  difference  equation  becomes 

(3M+2)TM  "  2TN  =  (3M~2)TM  +  2TN  (l6) 

4.  Method  of  solution 

Because  of  the  nonlinearities ,  an  iterative  solution  is  used. 

However,  since  the  coefficient  matrix  is  tridiagonal,  it  is  convenient 

and  expeditious  to  use  a  pivotal  method  several  times  rather  than  to 
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use  a  purely  iterative  method  such  as  that  of  Gauss-Seidel.     We  nay 
think  of  the  equations  in  the  form  (cf.   Reference  l'l  ) 


a 


iTi-l  +  biTi  +  ciTi+l  =  di  Ci-l,2....fNtM)   (17) 


with  a,  -  c..  =  0.  Then  we  form 

o*  =  0;  o*  =  V(bi_al0I-l5   (1=2,. ..,!!)  (18a) 

d*  =  0;  «  =  (drtd!.1)/(bi-aic!.1)  (1=1>2 M)    (lflb) 

and  obtain  the  advanced  temperatures  as 

TM  =  dM;     Ti  =  di  "  cFi+l       (i=N,N-l,...,2.1)  (19) 

Because  of  the  appearance  of  the  advanced  temperatures  in  the  coef- 
ficients, several  iterations  are  required. 

5.  Digital  computer  program. 

The  digital  computer  program  discussed  in  appendix  C  hereto 
implements  the  theory  developed  above.  The  main  program  is  called 
PIPETEM.  The  properties  of  the  fluid  are  determined  by  subroutines 
CSUBP,  PRAND,  and  ABSMU.  These  are  called  by  subroutine  AITCH  which 
determines  the  surface  heat  transfer  coefficient  h.  "The  only  fluid 
whose  properties  are  incorporated  in  the  present  program  is  superheated 
steam.  The  properties  of  the  metallic  pining  material  are  determined 
by  subroutines  METAL,  COEXP,  and  EYUNG.  All  these  subroutines  will 
be  discussed  further  in  appendix  A  hereto. 

6.  Properties  of  the  solution. 

Having  determined  the  advanced  temperatures  T. ,  (i=l,. . . ,M) , 
certain  significant  properties  of  this  solution  are  determined  in  a 
subroutine  called  AVERJ.  These  properties  are  TAV,  DTI,  DT2,  SL, 
and  SR  which  will  be  defined  and  described  at  this  point. 
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TAV,  DTI,  and  DT2  are  rigorous  counterparts  to  "average  temper- 
ature," "delta-tee-one,"  and  "delta-tee-two"  given  in  various  piping 
and  vessel  codes  for  the  estimation  of  thermal  stress.  As  will  be 
pointed  out  in  appendix  B,  for  uniform  cylinders  with  axial  symmetry 
of  loading  and  temperature,  the  proper  stresses  to  be  used  in  code 
evaluations  may  be  detsrmined  directly  without  the  use  of  these  quan- 
tities; however,  their  use  is  widespread  and  they  are  calculated  for 
comparison  purposes.  The  present  writer  is  responsible  for  the  defi- 
nitions given  in  the  codes.  For  simplicity,  these  definitions  were 
based  on  flat  plate  geometry  and  this  is  quite  reasonable  for  pipes 
except  those  having  unusually  thick  walls.  That  is,  in  these  defin- 
itions, pipe  diameter  was  assumed  to  be  quite  large  comparted  to  pipe 
wall  thickness.  In  subroutine  AVERT  this  assumption  is  not  made.  The 
appropriate  definitions,  based  on  truly  cylindrical  geometry,  are  given 
below,  followed  by  demonstrations  that  they  reduce  to  the  code  defin- 
itions as  pipe  radius  goes  to  infinity  while  maintaining  constant  wall 
thickness . 

TAV,  also  designated  as  T~,  is  the  averaps  metal  temperature  de- 


(20) 


where  a  =  inside  radius,  b  =  outside  radius,  h  =  b-a  =  wall  thickness 

(do  not  confuse  with  earlier  use  of  h  as  heat  transfer  coefficient), 

and  c  =  (b+a)/2  =  mean  radius.     If  we  define  y  =  r-c,  then 
r+h/2  r+h/2 

T  =  (lAi)\         (l+y/c)Tdy  -*  (1/h)  Tdy  as  c  -  °°  (21) 

J-h/2  J-h/2 

and  this  is  the  definition  given  in  the  codes. 
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fined  as 

b 

rb 

fb 

f  = 

*  /TdA//dA  = 

2  it 

Trdr/27T 
a 

rdr  = 

a 

(lAic) 

Trdr 
a 

We  next  define  a  quantity  which  can  be  termed  the  thermal 
moment,  viz. : 


U  =     r(T-T)dA  =  2tt 


r2(T-T)dr  =  2ttc 


+h/2 

(l+y/c)2(T-T)dy       (22) 


a  J -h/2 

Using  the  definition  of  T,  this  can  be  simplified  to  the  form 
r+h/2 


M  =  2ttc 


(l+y/c)Tydy  -  TTTh3/6  (23) 

-h/2 


The  next  step  in  this  development  is  to  consider  the  average  and 
the  thermal  moment  of  a  linearly  varying  temperature  distribution 
which  has  variation  V  from  r  =  a  to  r  =  b.  Such  a  temperature  is 
given  by  T  =  Vr/h  for  which  we  find 

rb 

?=  (V/h2c)  r2dr  =  (V/3h2c)(b3-a3)  (2*0 

b 
M  =  27t[    r2(Vr/n-T)dr  =   (irVch2/6)(l-h2/12c2)  (2r;) 

-"a 
It  requires  considerable  manipulation  to  obtain  the  last  form. 

By  equating  these  two  values  of  M,  one  evaluates  the  (equiva- 
lent)  linear  variation 

r+h/2 
V  =  [12  (l+y/c)Tydy  -  Th3/c]/h2(l-h2/12c2 ) 

J~h/2  ,+h/2 

-  (12/h2)  Tydy         (26) 

*  -h/2 

as  c  -*■  °°.     The  limiting  values,  as  c  ■*  °°,     are,  in  each  case,  the 

definitions  given  by  the  various  piping  and  vessel  codes.     The  present 
subroutine,  AVERJ,  employs  the  more  general  definitions  based  on 
cylindrical  geometry  with  finite  radii.     The  variation  V  is  also 
denotes  as  ATj   in  the  codes  and  as  DTI  in  the  output  of  the  present 
program. 
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In  the  codes  the  quantity  AT2  is  defined  as  the  nonlinear 
residual  temperature  at  the  inner  or  the  outer  radius  (whichever 
gives  the  greater  result)  after  the  constant  part,  T,  and  the  linearly 
varying  part,  V,  are  subtracted.  In  the  outnut  of  the  present  program 
this  quantity  is  denoted  as  DT2  and  defined  as 

DT2  =  I4ax{  |  T^T-OTl/2 1  ,  |  T  -TfDTl/2 1  }  ( 27 ) 

The  quantities  ATX  and  AT2  are  intended,  in  the  codes,  to  permit 
estimating  significant  thermal  stresses  in  general  pressure  vessel 
components.  However,  in  the  particular  case  of  uniform  hollow  cylin- 
ders, the  stresses  resulting  from  radial  temperature  gradient  are 
definitely  classified  as  "local  thermal  stresses"  and  there  is  no 
actual  need  to  separate  the  actual  temperature  distribution  as  has 
been  indicated  above.  Instead,  the  thermal  stresses  can  be  calculated 
directly,  and  this  is  done  in  subroutine  AVERT.  The  quantities  SL 
and  SR  are  the  values  of  circumferential  (or  axial  —  they  are  equal) 
stress  due  to  the  temperature  distribution  at  inside  and  outside  metal 
surfaces,  respectively,  under  the  assumption  of  zero  net  axial  force 
and  bending  and  twisting  moment.  The  latter  contributions  may  be  de- 
termined separately  by  use  of  a  conventional  piping  flexibility  analysis, 

Because  of  the  variability  (with  temperature)  of  the  Young's 
modulus  and  of  the  coefficient  of  thermal  expansion,  the  calculations 
are  not  trivial.  Discussion  of  this  matter  is  given  in  appendix  B 
hereto. 
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APPENDIX  A       Discussion  of  the  subroutines 

Subroutines  CSUBP,  PRAND,  and  ABSMU  receive  as      irbuts  the  tem- 
perature  (°F)  and  the  pressure  (psia)   of  the  steam  flowing  in  the  pine 
and  produce  as  outputs  the  values  o£  c     (specific  heat  at  constant  pres- 
sure), Up  (Prandtl's  number) ,  and  u  (absolute  viscosity)   o^  the  steam. 
The  calculations  are  by  polynomials  in  T  and  P  the  coefficients  of  which 
were  obtained  by  least  squares  analysis  of  data  points  read  from  the 
graphs  on  pages  293,  29^,  and  297  of  the  ASP-IE  Steam  Tables  (      ).     It  is 
implicitely  assumed  that  the  steam  is  saturated  or  superheated.     If  the 
steam  is  wet,  the  proper  values  are  orobably  not  obtained. 

Sixty  number  pairs  were  used  in  dtermining  the  coefficients  in 
CSUBP,  fifty-four  pairs  for  PRAND,  and  twenty- four  pairs  for  ABSMU.     The 
only  type  of  accuracy  evaluation  performed  v/as  to  observe  that  each  o^ 
the  subroutines  was  able  to  generate  the  input  number  pairs  within  one 
or  two  percent  (generally  much  closer  than  this).     Such  accuracy  is 
greater  than  that  implicit  in  convective  heat  tranfer  theory  and  certainly 
at  least  as  great  as  the  accuracy  with  which  the  input  values  were  read 
from  the  graphs. 

Subroutine  AITCH  determines  the  film  heat  transfer  coefficient  h 
by  use  of  the  Colbum  formula  (cf.  Giedt  (  8  )) 

h  =  NgjQc  (A-l) 

where  G  is  mass  per  unit  time  per  unit  area,  i.e., 

G  =  ^Aflow  (A-2) 

and  NLm  is  the  Stanton  number, 
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N_  =  0.023  iC2/3  rC1^  (M) 


ST       "'"'-'  "PR      "re 
In  this  formila  NpR  is  the  Prandtl  number  and  IJRF  is  the  Reynolds  number 

NRE  =  Gd/U  ^'^ 

where  d  is  the  pipe  inside  diameter.     The  Colburn  formula  is  one  which 

is  recommended  for  cases  of  well  developed  turbulent  flow     at  fairly 

hipft  Reynolds  numbers.     It  seems  to  be  an  appropriate  choice  except 

perhaps  for  cases  in  which  the  operating  transient  involves  a  very 

great  decrease  in  mass  flow  rate.     In  such  cases    Q indeed  approaches  zero 

but  not  necessarily  in  the  way  represented  by  the  Colburn  formula. 

Incidentally,  subroutine  AITCH  avoids  a  computational  difficulty  as 

G  goes  to  zero  by  combininr,  the  equations  in  the  form 

h  =  0.023  cp  G0*8  (u/d)0,2  /  N^3  (A-5) 

The  specific  heat  c  is  evaluated  at  the  fluid  bulk  temperature 
but  NpR  and  p  are  evaluated  at  the  "film"  temperature  which  is  the 
arithmetic  average  of  fluid  bulk  temperature  and  inside  wall  metal 
temperature.  Since  the  latter  is  an  unknown  in  the  calculations,  an 
iterative  evaluation  must  be  used. 

Subroutines  METAL,  COEXP,  and  EYUHG  obtain  pertinent  properties 
of  the  metal  oipe  material  corresponding  to  input  values  of  the  metal 
temperature.  These  subroutines  presently  provide  data  for  only  one 
particular  material,  the  standard  low  chrome-mo ly  high  temperature 
alloy  designated  by  the  number  22  —  pipe  is  designated  as  P22.  Sub- 
routine METAL  provides  values  of  thermal  conductivity  (k),  the  rate 
of  variation  (k'  =  dk/dT) ,  and  the  thermal  diffusivitv  using  poly- 
nomial representations  of  material  test  data. 
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Although  Table  1-4.0  of  the  ASME  Nuclear  Components  Code  (2) 
suggests  using  the  same  values  for  thermal  conductivity  and  thermal 
diffusivity  for  P22  as  for  low  carbon  steel,  data  given  in  the  T^RC 
Data  series,  References  (15)  and  (16),  indicate  a  sufficiently  great 
dependence  on  alloying  constituents  that  we  have  preferred  to  use  the 
latter.  Thus  we  get  substantially  lower  values  of  conductivity  and 
diffusivity  than  the  ASME  tabulation  gives.  The  actual  data  employed 
in  constructing  subroutine  METAL  are  the  data  points  from  curve  19 
page  1158  of  Ref.  (15)  and  curve  2  page  3^3  of  Ref.  (16).  It  should 
be  noted,  however,  that  it  would  be  a  very  simple  matter  to  modify 
the  subroutine  METAL  to  use  the  ASME  data  rather  than  the  TPRC  data. 

Subroutines  COEXP  and  EYUNG  are  used  (by  subroutine  AVERG)  to 
provide  the  coefficient  of  thermal  expansion  and  the  Young's  modulus 
corresponding  to  the  given  input  temperature.  The  data  is  from 
Appendices  B  and  C  of  the  Power  Piping  Code,  Ref.  (3),  and  linear 
interpolation  is  used.  Poisson's  ratio  is  assumed  to  have  the 
constant  value  0.3. 

All  subroutines  as  well  as  the  main  program  use  Elnglish  custom- 
ary units,  specifically  inches  (rather  than  feet)  and  seconds  (rather 
than  minutes  or  hours ) . 
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APPENDIX  B    Stress  Calculations 

Assuming  what  is  called  a  "quasi-static"  situation,  i.e., 
neglecting  forces  associated  with  accelerating  the  mass  of  the  oipe 
as  its  thermal  expansion  causes  displacement  of  its  material  particles, 
and  assuming  also  that  the  difference  between  adiabatic  and  isothermal 
material  constants  is  of  neglegible  significance,  we  nevertheless  are 
faced  with  a  nontrivial  problem  in  determining  the  state  of  stress 
which  results  from  the  temperature  changes  discussed  in  the  body  of 
this  report. 

Two  material  constants,  E  (Young's  modulus)  and  ft  (the  coefficient 
of  linear  expansion  —  do  not  confuse  with  thermal  diffusivity)  are  of 
significance  in  connection  with  stress  evaluation.  Roth  vary  with 
temperature.  Further,  as  will  be  discussed  below,  a   varies  with  stress 
state.  A  number  of  writers  have  dealt  with  stresses  in  tubes  due  to 
changes  in  temperature,  among  them  Hilton  (  9)  and  Chang  and  Chu  (7). 
The  analysis,  in  effect,  reduces  to  the  following. 

ap  =  [E/r2(l-v)][(r2-a2)  J(b)/(b2-a2)  -  J(r)]  (B-l) 

aQ  =  [E/r2(l-v)][(r2+a2)  J(b)/(b2-a2)  +  J(r)  -  r2aT(r)]    (B-2) 
6 

a    =  Eg  +  [E/(l-v)][2vJ(b)/(b2-a2)  -  tfT(r)]  (B-3) 

z    z 


where  „ 

J(r)  = 


a 
Presuming  that  the  axial  force  P  is  zero,  we  get 

e  =  2  J(b)/(b2-a2)  (B-5) 

z 

and  can  thus  evaluate  a   .  For  all  temperature  distributions  which 
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result  from  external  heating  (i.e.,  from  heat  applied  to  the  inner  or 
outer  surfaces  but  not  generated  within  the  wall  itself)  the  extremal 
stresses  will  be  at  r  =  a  or  at  r  =  b.     (With  external  insulation,  it 
is  plausible  that  the  extreme  must  be  at  r  =  a. )       V/e  thus  have  a     =  0 
at  r  =  a  and  at     r  =  b,  and  we  also  have 

a0  =  oz  =  [E/(l-v)][2J(b)/(b2-a2)  -  nff(a)]     at  r  =  a         (B-6) 
aQ  =  oz=  [E/(l-v)][2J(b)/(b2-a2)  -  aT(b)]     at  r  =  b         (B-7) 

In  computing  J(r)  and  J(b),  the  a  under  the  sign  of  integration 
is  the  "instantaneous"  value  appropriate  to  T(r).     It  is   (nresumably) 
the  "zero-stress"  value  of  a  which  is  tabulated.     Obviously,  however, 
during  the  temperature  changes  the  various  elements  of  the  pine  wall 
are  not  stress-free     and  the  question  arises  of  how  it  is  that  we 
should  use  the  zero-stress  value  of  a.     The  explanation  is  surely  not 
well  known  (this  question  is  not  even  considered  in  any  of  the  standard 
references  on  thermal  stress  analysis);  the  explanation  which  follows 
was  given  by  the  writer  in  a  report,  dated  November  i960,  entitled 
"Reflections  on  Piping  Flexibility  Analysis,"  which  was  intended  to  ex-' 
plain  the  basis  of  procedures  incorporated  in  the  "MDC  Document"  of 
the  Mechanical  Design  Committee  of  the  ANSI  (then  ASA)  B31  Piping  Code 
Committee. 

This  explanation  considers  a  one  dimensional  constrained  thermal 
expansion  problem  of  a  carbon  steel  bar  of  length  L  which  at  70  F 
fits  precisely  between  two  rigid  and  immovable  anchors;  that  is,  at  70  F 
there  is  no  stress  and  no  gap.     The  temperature  is  increased  to  300  v. 
One  is  to  determine  the  resulting  compression  stress,  nresuming  that 
there  is  no  tendency  to  buckle  laterally.     As  fundamental  we  take  the 
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relation 

de  =  da/E  +  adT  (U-8) 

and,  since,  in  this  problem  de  =  0,  we  have 

da  =  -EmdT  (B-9) 

r300 

E(T)a(T)dT  (B-10) 

70 


aP=  " 


where  subscript  (p)  indicates  the  "final"  condition. 

Appendices  B  and  C  of  MSI  B3 1.1-1973,  Ref.  (   ),  give  the 
following  data  for  carbon  steel  pipe: 

T      E'10"6      rx-106       En 


70  27.  6.07  169.35 

200  27.7  6.38  176.73 

300  27.4  6.60  180.84 

Using  trapezoidal  integration,  we  evaluate  equation  (B-10)   to  get 

ap  =  -(169. 35+176. 73)(130/2)+(176. 73+180. 84)(100/2) 

=  -4037^  psi  (B-ll) 

(This  result  is  incorrect  as  we  will  show.) 

An  alternate  calculation  obtains  the  answer  by  considering  one  anchor 

to  be  removed  and  letting  the  bar  expand  freely  («r»-stress  condition) 

from  70  F  to  300  F  and  then  supplying  sufficient  compressive  load  to 

restore  the  bar  to  its  original  length.  This  calculation  gives 

(T)dT  =  [(6.07+6.38)(65)+(6.38+6.60)(50)]-10-t) 


a 
e  = 


70 


=  1.4583-10-3  (B-12) 

a„  =  -(27.4»106)(1.4583-10~3)  =  -39956  psi  (B-13) 

r 

which  is  a  slightly,  but  definitely,  smaller  value. 

Now  it  would  indeed  be  disturbing  if,  without  any  inelastic 
behavior,  the  final  stress  state  were  to  denend  on   the  sequence  of 


-15- 


load- temperature  application.  This  would  indeed  be  a  new  type  of 
thermal  ratcheting.  However,  we  will  now  show  that  the  second  cal- 
culation gives  the  correct  result  for  any  load- temperature  sequence, 
based  upon  the  postulate  that  the  final  stress  state  is  independednt 
of  sequence  and  depends  only  on  material  pronerties ,  the  final  tem- 
perature state,  and  the  final  strain.  In  other  words,  we  nresume  that 
equation  (B-  8)  represents  an  exact  differential  so  that 

ep  =  f(ap,Tp) 
If  this  is  the  case,  we  have 

(B-14) 

(B-15a,b) 

(B-16) 

(we  have  written  d  rather  than  9  on  the  right  since,  by  its  definition, 
E  does  not  depend  on  a . ) 

Thus,  if  E  varies  with  T,  then  a  must  vary  with  a;  i.e.,  a  = 
ot(T,a).  The  tabulated  values  we  use  are  values  of  a(T,0)  and  these 
were  appropriate  in  the  second  numerical  calculation  but  not  in  the 
first.  Thus,  for  our  original  problem  we  should  integrate 

do  =  -E(T)rx(T,a)dT  (B-17) 

taking  into  account  the  dependence  of  a  upon  a.  We  have 

a(T,a)  =  a(T,0)  +  J  (|g)d0  =  a(T,0)  +  a^(|)  (B-lB) 

so  that 

da  =  -E(T)[a(T,0)  +  (a/E2)(dE/dT)]dT 
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9   (de. 

9  /  de,. 
9cT9T; 

and,  since 

9e 

3T  -  a» 

9e  .  1 
9a       E 

we  have 

9a       d  /l 
9a  ~  dT^E 

) 

=  -E(T)oc(T,0)dT  -  adE/E  (3-19) 

Thus, 

-a(T,0)dT  =  (Eda  -  cdE)/E2  =  d(a/E)  (13-20) 

Performing  this  integration,  we  get 

r300 
(a  /E  )  =  -         a(T,0)dT  =  -1.4583-10"^ 
HO 

ap  =  (-U»583'10"3)(27J4-106)  =  -39956  psi 
as  in  the  second  calculation. 

In  the  more  complicated,  multi-axial  stress  situation  with  which 
this  monograph  is  primarily  concerned,  the  argument  takes  the  following 
form.  We  imagine  any  element  (dr,rd  ,dz)  of  the  pipe  to  be  separated 
from  the  adjacent  elements  and  to  experience  stress- free  thermal  ex- 
pansion resulting  in  extensional  strains  of  the  form 

e  =    a(T,0)dT  (B21) 

Then  surface  forces  are  applied  to  the  element  to  cause  it  to  assume 
a  form  such  that  all  elements  can  be  fitted  to  gether  in  the  final 
stressed  state.  Since  equations  (B-l)  through  (B-7)  were  derived  by 
satisfying  the  equilibrium  and  constitutive  laws,  they  indeed  describe 
the  final  state. 

One  final  remark  needs  to  be  made  concerning  stress  calculations 
for  code  purposes.  Although  the  so  called  "simplified"  methods  of  an- 
alysis provide  for  stress  components  corresponding  to  the  temperature- 
like quantities  ATt  and  AT2 ,  it  is  very  important  to  note  that  in  the 
specific  case  of  a  uniform  pipe  or  cylindrical  vessel  the  stress  state 
resulting  from  a  radial  temperature  gradient  does  not  have  to  be  separ- 
ated into  "secondary"  and  "local"  parts  but,  instead,  the  actual  stress 
may  all  be  considered  as  "local."  The  word  "local"  is  quite  important 
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in  considering  allowable  states  of  stress.  A  definitive  statement  to 
tnis  effect  is  to  be  found  on  the  fifth  page  of  the  ASME  "Criteria" 
document,  Ref.  (  4 ),  viz: 

"A  special  exception  to  these  general  rules  is  the  case  of  the 
stress  due  to  a  radial  temperature  gradient  in  a  cylindrical 
shell.  It  is  specifically  stated  in  II-'ll2(m)(2)(6)  of  Section 
III,  and  in  4-112(1) (2) (6)  of  Appendix  'I  of  Division  2,  that 
this  stress  may  be  considered  a  local  thermal  stress.  In 
reality  the  linear  portion  of  this  gradient  can  cause  defor- 
mation, but  it  was  the  opinion  of  the  Special  Committee  that 
this  exception  could  be  safely  made." 
Thus  the  stresses  calculated  according  to  equations  (B-6)  and  (B-7) 
of  this  appendix  themselves  directly  give  the  significant  thermal 
stresses  due  to  the  radial  gradient,  and  these  are  indeed  local 
thermal  stresses. 
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APPENDIX  C  Instructions,  Data  Preparation,  Typical  Problem 

Instructions  for  the  use  of  the  programs  and  preparation  of  in- 
put data  are  given,  in  the  form  of  comments,  in  the  initial  portion 
of  the  program  itself.  This  portion  of  the  program  reads  as  follows: 
Temperature  distribution,  etc.,  in  uniform  pipe  with  exterior 
insulation,  convectively  heated  by  fluid  contents.  The  basic 
procedure  is  a  finte  difference  treatment  of  the  one  (space) 
dimensional  partial  differential  equation  in  cylindrical  co- 
ordinates. Metal  properties  are  assumed  to  vary  with  temper- 
ature. Surface  heat  transfer  coefficient  is  determined  by 
Colbum's  relationship.  The  fluid  content  In  this  version  is 
assumed  to  be  superheated  steam  and  the  pipe  material  is  P22. 

Input  is  by  means  of  two  kinds  of  cards.  The  first  card  reads 
OD,  WT,  TINIT,  N,  MIT,  and  NPROB  according  to  FORMAT(3FiO. 0,315) . 
OD  is  the  outside  diameter  in  inches.  WT  is  the  wall  thickness 
in  inches.  TINIT  is  the  uniform  initial  metal  temperature  in 
degrees  Fahrenheit.  (For  nonuniform  initial  temperature,  see 
below.)  N  is  the  (even)  number  of  equal  subregions  into  which 
the  wall  thickness  is  divided.   (If  input  N  is  odd,  it  is  round- 
ed up  to  the  next  even  number. )  NIT  is  the  number  of  iterations 
employed  each  time  step  in  order  to  solve  the  nonlinear  matrix 
equation.  NPROB  is  an  integer  identification  number. 

If  TINIT  <  -500.,  the  program  looks  for  nonuniform  initial  metal 
temperature,  reading  enough  cards  according  to  t?ormat(  7F10 . 0 ) 
to  provide  the  (N+l)  required  initial  temperatures. 
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Subsequent  cards  read  I,  NTEES,  TIITL,  TIM2,PHI1,  PHI2,  PRES1, 
PRES2,  FL0W1,  and  FL0W2  according  to  FORMAT(2I5,8F8.0) .     I  is 
an  identifying  serial  nurrber  which  is  useful  for  ordering  the 
data  but  which  is  not  used  in  the  nrogram.     NTEES  is  the  nurrber 
of  time  sub intervals  into  which  the  large  time  interval  from 
time  =  TIM1  until  time  =  TIM2  is  to  be  divided.     PHI1  and  PHI2 
are  values  of  the  fluid  temperature  at  times  TIM1  and  TTM2 
respectively,  This  temperature  is  assumed  to  vary  linearly 
with  time.     Similarly  for  PRES  (fluid  pressure  in  psig)  and 
FUM  (fluid  flow  rate  in  pounds  per  second). 

The  data  deck  for  a  problem  terminates  with  a  card  containing 
10000  in  the  first  five  positions.     A  number  of  problems,  each 
having  a  data  deck  as  described  above,  may  be  stacked.     The 
sequence  of  problems  is  normally  terminated  by  a  card  contain- 
ing -1.   in  the  first  three  positions. 

Some  thought  was  given  to  providing  an  option  which  would  permit 
several  time  intervals  to  be  calculated  per  output  line  of  print. 
It  was  decided  that  not  only  would  this  be  a  needless  complication 
but  also  that  frequent  output  is  appropriate  and  useful  during  those 
periods  when,  for  one  reason  or  another,  small  time  increments  are 
employed. 

It  is  easily  possible  to  make  minor    internal  changes  so  that  the 
nodal  temperatures,  rather  than  the  results  calculated  in  AVER J, 
are  printed.     This  is  useful  when  comparing  results  with  those  ob- 
tained by  other  methods  of  calculation. 
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The  program  listing  is  not  given  here  for  reasons  of  space  economy. 
Those  who  are  seriously  interested  in  installing  the  nrogram  should 
communicate  directly  with  the  writer  to  make  appropriate  arrangements. 
The  program  is  written  in  FORTRAN  using  none  other  than  quite  primi- 
tive capabilities.     It  has  been  run  on  an  IBM  360/67  requirinr  only 
56K  core  for  execution  (both  G  and  II  compilations).     For  a  stack  or 
thirty  two  problems,  involving  a  total  of  5240  time  steps,  with  N  = 
10   (10  layers)   and  NIT  =  5  (five  iterations),  time  requirements  were 
as   follows :  G  Compiler        H  Compiler 


Compile 
Link-edit 
Execute 
Total 


14.38s  23.20s 

1.87s  1,74s 

I4ml9.48s  llml7.03s 

I4m35.73s  llm4 1.97s 
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APPENDIX  D  Program  verification 

By  artificially  causing  the  program  PIPETEM  to  operate  with  constant 
values  of  diffusivity,  conductivity,  and  surface  heat  transfer  coefficient, 
it  may  be  used  to  solve  problems  for  which  "classical"  solutions  are  known. 
The  program  was  initially  verified  by  taking  O.D.  =  10000.,  W.T.  =  1.,  and 
comparing  results  with  data  points  taken  from  graphs  given  by  McNeill  and 
Brock  (12),  When  a  number  of  such  comparisons  checked  within  the  accuracy 
with  which  the  graphs  could  be  read,  a  more  exacting  test  was  made.  This 
is  described  in  what  follows. 

Ozisik  (13)  gives  solutions  involving  true  cylindrical  geometry,  ^or  a 
uniform  cylinder  with  inside  radius  a  and  outside  radius  b,  insulated  on 
its  outer  surface,  and  conveetively  exposed,  with  constant  surface  heat 
transfer  coefficient  h,  to  an  internal  fluid  having  "ramp"  temperature 
<J>  =  At  (where  A  is  a  constant  and  t  denotes  time) ,  and  having  initial 
temperature  distribution  T(r,0)  =  0  at  time  t  =  0,  the  solution  may  be 
written  as 


Xr,t)  =  7  t  (t)ib  (p)/A  \j)   (1) 
'      S  m   rm    mm 
m=l 


where 


t  (t)  =  2ABrt[exp(-Y  t)  -  1  +  y  t]/a2 
m  m  'm 

Y     =  nt(u  /a)2 
m  m 

ib  (p)  =  Yn(py  )/Y,(py  )  -  J0(pu  )/J,(py  ) 
m  °       m       i       m  °       m       l       m 

A     =  y  2{[p^  (p)/ib  (l)]2   -  1  -(B/p   )2} 
m        m         rm         rm  m 

B  =  ahA  (Biot's  number),   p  =  r/a,  o  =  b/a 

and  the  numbers  u     are  the  zeros  of  the  function 
m 

[BJ0(v0  +  yJjCyJjYjCpu)  -  [BYQ(p)  +  yYl(y)]JI(py) 

The  symbols  J.   and  Y.   denote  Bessel  functions. 
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It  proved  to  be  quite  difficult  to  devise  a  prolan  to  nerform 
these  calculations.  An  apparent ly  different  solution,  given  by  Luikov 
(lO),  was  dealt  with  first,  but  debusing  was  not  completed  successfully. 
A  later  work  by  Luikov  (11),  in  Russian,  gave  the  form  v/e  actually  used, 
citing  Ozisik  (13)  as  source.  Both  single  and  double  precision  versions 
of  an  implementing  program  were  constructed  and  debugged.  The  double 
precision  version,  much  slower  than  its  counterpart,  partly  because  of 
the  use  of  a  "home-made"  DP  Dessel  function  subroutine,  was  used  to 
verify  that  the  single  precision  version  provided  adequate  accuracy. 
In  all  such  checks,  results  from  the  two  propxams  agreed  to  five  sig- 
nificant digits. 

The  test  problem  employed  to  check  the  integrity  and  accuracy 
of  PIPETEM  was  as  follows:  a  =  4  inches,  b  =  5  inches,  A  =  1. ,  B  =  1. , 
a  =  1. ,  k  =  1. ,  k'  =  0. ,  h  =  0.25.  Analytical  calculations  were  made 
for  t  =  0,  1.,  2.,  3.,  ^.,  and  5.  PIPETEM  calculations  were  made  for 
t  =  o.  to  t  =  5.  at  intervals  of  0.1.  In  PIPETEM  we  took  N  =  10  and 
HIT  =  5.  (V/e  wasted  a  few  seconds  by  not  usinf  NIT  =  1.  This  problem 
is  completely  linear  and  no  iteration  is  necessary.)  In  the  analytic 
solution  initially  we  took  account  of  ten  terms  in  the  summation. 

The  results  were  so  close  that  differences  can  not  be  displayed 
graphically.  However,  the  differences  were  sufficiently  great  as  to 
demand  an  explanation.  Truncation  error  in  the  analytic  solution  was 
suspected.  Accordingly,  the  analtic  solution  was  repeated  three  more 
times  using  five,  twenty,  and  forty  terms.   (The  forty  term  solution 
was  very  extravagant  of  computer  time.)  Plotting  typical  analytic 
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results  versus  reciprocal  number  of  terms  clearly  indicated  that  indeed 
the  analytical  solution  did  suffer  from  truncation  error  and  that  extra- 
polation did  indeed  lead  to  the  result  given  by  PIPETEF4.     Since  the  theo- 
retically "exact"  solution  has  by  far  the  greater  error,  it  is  not  easy 
to  estimate  the  accuracy  of  the  PIPETEh  results,  but  in  this  problem  the 
error  appears  to  be  only  Htt.lt  n<*rothan  one  part  in  one  thousand.  The 
following  tabulation  gives  calculated  results  for  T(0,5),  for  example. 

N  =  5         N  =  10         N  =  20         N  =  40         N  =  00        PIPETEM 
2.05694     2.08627       2.09991       2.10647       2.11496       2.11220 
Table  D-l     Calculated  results  for  T(0,5) 

The  first  four  tabulated  results  were  obtained  by  the  single  precision 
analytic  solution.   (The  double  precision  version  also  gave  2.0R627  for 
N  =10.)     The  extrapolated  value,  for  N  =  00  ,  was  obtained  by  assuming 
that  T  =  A  +  Bx  +  Cx2  +  Dx3,  where  x  =  5/II,  leading  to  the  evaluations: 
A  =  2.114958096,  D  =  -.073526667,  C  =  .049093333,  D  =  -.033584762.     The 
difference  between  the  extrapolated  analytic  value  and  the  value  provided 
by  PIPETEM  is  0.00275  which  is  0.13$  of  the  average  of  these  values. 

As  an  ultimate  check  on  the  validity  of  the  program,  an  artificial 
nonlinear  problem  was  made  up  and  solved  numerically.     The  conditions 
are:   inside  radius  a  =  1,  outside  radius  b  =  2,  conductivity  k  =  T  (so 
that  k'   =1),  heat  transfer  coefficient  h  =  15e  ,   fluid  temperature  <t> 
=  6e  ,  initial  temperature  distribution  T(r,0)  =  r  +  4/r.     The  diffus- 

ivity  a  is  a  function  of  both  temperature  and  time,  viz: 

-t 


a 


=  Y/[l+(4-Y)2/(4+Y)2],  Y  =  (Z-/Z2-l6)2/4,   Z  =  Te 


It  is  easy  to  verify  that  T(r,t)  =  e  (r  +  4/r)  is  a  solution  of  this 
problem.  Although  matters  of  uniqueness  of  solution  to  nonlinear  sys- 
tems are  difficult  to  deal  with,  it  happens  that  PIPETEM  indeed  yields 
the  solution  cited.        _2fi 


Attention  was  focussed  on  the  values  T  =  T(l,t*)  =  100.000 
and  TB  =  T(2,t*)  =  80.000,  where  t*  =  log  20.  Percent  errors  are 
shown  in  Table  D-2  corresponding  to  the  128  evaluations  indicated, 
values  of  r  being  1,  2;  values  of  N  being  2,  4,  6,  10;  values  of 
NIT  being  1,  2,  5,  10;  and  values  of  CTEE  being  t*/5,  t*/20,  t*/50, 
t*/200. 


DfEE 


t*/5 


NIT 


N  =  2 


N  =  4 


II  =  6 


M  =  10 


t*/20 


t*/50 


t*/200 


1  -21.5V-l4.04  -14.70/-7.173  -9.532/-2.487  -3.372/+3.679 

2  -2.726/+1.479  +.3159/+3.258  +.0102/-. 0009  -. 7222/-2. 6^6 
5  -2.592/-6.721  -1.78V-6.331  -.9409/-3.697  +1.2R9/+3.447 

10  -.8019/+4.055  +.5816/+3.232  -.0216/+. 3814  -2.023/-6.768 

1  -2.042/+2.726  +1.085/+4.366  +1.840/+4.425  +2.569/+4.810 

2  -1.962/-. 9932  -1.070/-1.907  -1.010/-2.239  -l.inV-2.59^ 
5  -1. 648/+. 8822  -.3843/+. 7397  -.IO63/+.5389  +.0671/+.5541 

10  -1.735/-. 1605  -.5519/-. 0690  -.2782/-. 0115  -.1289/-. 1902 

1  -1.379/+2.520  +.2103/+2.209  +.6620/+2.l8I5  +.9698/+2.209 

2  -1.840/-. 1536  -.7005/-. 5568  -.4543/-. 5960  -.3513/-. 6870 
5  -1.714/+.4576  -.5162/+. 1779  -.2334/+. 1084  -.0781/+. 0837 

10  -1. 723/+. 3710  -.5285/+. 0778  -.2466/-. 0069  -.0916/-.0043 

1  -1.631/+1.020  -.3372/+. 7247  -.0229/+. 6389  +.1658/+.6182 

2  -1. 734/+. 3131  -.5426/+. 0234  -.2618/-. 0232  -.1089/-. 0567 
5  -1. 725/+. 3715  -.5287/+. 0857  -.2464/+. 0365  -.0910/+. 0131 

10  -1. 725/+. 3713  -.5288/+. 0489  -.2465/+. 0351  -.0911/+.0114 


Table  D-2     Percent  errors  in  evaluation  of  T.  and  TR  in  test  problem. 
Error  in  T.   (Tg)  appears  to  left  (right)  of  slash  (/). 

It  may  be  noted  that  the  accuracy  is  quite  sufficient  for  engineering 
purposes  if  the  number  of  layers  N  =  4,  the  number  of  iterations  MIT 
=5,  and  the  total  time  interval  is  divided  into  20  parts.     It  may 
also  be  remarked  that  the  total  time  to  evaluate  all  128  results  was 
only  2  minutes  20  seconds. 
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APPENDIX  E   Heat  Transfer  at  Outer  Radius 

Following  the  completion  of  the  preceding  portions  of  this  renort, 
an  apDlication  came  to  attention  which  called  for  the  alimentation  of 
PIPETEM  capabilities.  Ordinarily  high  temperature  nip inn  and  fitting 
are  lagged  (insulated).  However,  in  some  cases  it  might  be  of  inter- 
est  to  follow  the  temperature  history  of  unlapged  pipes  when  subjected 
to  internal  fluid  transient  temperatures.  Also,  PIPETEM  is  obviously 
applicable  to  insulated  flanges  if  one  discounts  axial  conduction  and 
the  additional  thermal  resistance  of  the  flange-bolt  interfaces. 
Down-transients  (i.e.,  decreases  in  the  temperature  of  fluid  contents) 
lead  not  only  to  local  thermal  stresses  but  also  to  discontinuity 
thermal  stresses  if  there  are  significant  changes  in  wall  thickness, 
since  the  average  temperature  of  the  thicker  part  lags  (in  time)  be- 
hind that  of  the  thinner  part.  This  is  particularly  true  of  flange- 
pipe  assemblies.  Leaving  the  flange  unlagged  (but  with  lagging  on 
the  pipe  itself)  causes  the  average  temperature  of  the  flange  to  be 
lower  than  that  of  the  pipe  since  the  former  can  dissipate  heat  by 
convection  and  radiation.  Then,  when  a  down- transient  occurs,  the 
resulting  thermal  dislocation  stress  will  be  less  than  what  it  would 
have  been  if  the  flange  had  been  lapged.  The  calculation  of  the 
temperatures  in  an  unlagged  flange  really  calls  for  finite  element 
treatment  since  axial  conduction  and  heat  loss  ^rom  bolts,  nuts,  and 
flange  face  is  obviously  of  importance.  However,  a  useful  estimate 
of  the  effect  can  be  gained  by  considering  only  radial  flow,  with, 
perhaps,  an  artificially  enhanced  value  of  the  effective  convective 
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external  heat  transfer  coefficient. 

Accordingly,   PIPETEM  has  been  modified  to  account  for  a  convective 
condition  at  the  outer  boundary,  viz. 

k^B  =  h(-^TM)  E_l 

rather     than  B  =  0.     The  nth  finite  difference  equation  becomes 
%  +2  +  a+[2  +  ±  -  o^+(T^+).]}T;  -  2^ 

=   C6M-2)TM  +  2TN  -  a(Tr,f.!;)[2  +  i  -  naC?^,)] 

+  aV(2  +  ^+  oJpV)  E-2 

This  equation,  which  takes  the  place  of  equation  16  and  which  reduces 

to  equation  16  when  h  (external)  =  0,  may  be  comnared  term  by  term  with 

equation  14.  In  equation  E-2,  quantities  have  the  meanings  ascribed 

to  them  on  page  5>  but  the  evaluations  are  at  the  outer  radius.  Note 

that  all  signs  are  the  same  (as  in  equation  1^0  exceot  for  the  —  terms. 

The  quantity  ^  is  the  external  ambient  temperature  ( for  convection) . 

External  heat  transfer  can  take  place  both  by  convection  and  by 

radiation,  so  that  the  coeffcient  h  in  equation  E-l  is  the  sum  of 

two  coefficients 

h  =  h  +  h  E-l 

c    r 

The  first,  h   ,   is  the  actual  convective  coefficient.     We  have  presumed 
'     c' 

natural  convection  from  a  cylinder  into  still  air,  using  the  relation 

h  =  0.2'l(eAi)0,2Cj  E-2* 

on  page  218  of  reference    8 »     Here 

0  =   |TH  -  ij/|  E-5 

and  d  is  the  diameter  of  the  cylinder,  i.e.,  the  outer  diameter  of  the 
pipe.  The  "radiative-convective"  coefficient  h  is  an  equivalent  con- 
vective coefficient  to  account  for  radiative  heat  transfer  and  is  de- 
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fined  by  the  equation 

Q  =  aeA[(TM+C)^-(^-C)',]  =  hpA(TB-^)  E-6 

in  which  a  is  the  St  e  fan-Bolt  zmann  constant  (0.171>J  xl0~8  Btu/hr  ft2  ^k)3 
z   is  the  emissivity  (we  take  0.8  for  bare  steel),  A  is  area  (which  cancels), 
ijj  is  the  temperature  of  the  surrounding  with  which  radiative  heat  trans- 
fer takes  place  (this  may  be  and  usually  is  different  from  the  convective 
ambient  temperature  ijj),  and  C  =  460  is  the  constant  required  to  convert 
Fahrenheit  temperatures  to  Rankine  (absolute)  temperatures. 

The  value  of  h  determined  from  equation  E-3  is  in  units  of  feet  and 
hours  rather  than  inches  and  seconds  and  it  must  be  multiplied  by  l/(720)2 
before  being  used  in  PIPETEM.  There  is  provision  in  PIPETEM  for  using 

the  relation 

h  =  f  h  +  f  h  E-7 

c  c    r  r 

where  f  and  f  are  "augmenting  factors"  which  may  be  used  to  account 
for  the  fact  that  the  geometry  is  not  exactly  that  of  an  infinite  cylinder 
so  that  a  "face"  or  protuberances  such  as  nuts  and  bolts  nay  convect 
and/or  radiate  more  heat  than  would  otherwise  be  accounted  for. 

The  convective  ambient  temperature  ^  and  the  radiative  ambient 
temperature  \p     could  be  made  to  be  functions  of  time,  just  as  with  the 
internal  fluid  temperature  d>,  but  this  has  not  been  done  in  the  nresent 
implementation  in  which  ty   and  ijJ  are  taken  to  be  constant. 

The  input  (and  the  input  instructions)  have  been  modified  slightly 
to  accommodate  this  added  facility.  For  each  problem  an  additional  card 
is  inserted  after  the  first  card  described  previously.  If  this  card  is 
blank,  perfect  insulation  is  assumed  at  the  outer  surface.  Otherwise, 

the  card  inputs  ^,  tyt   and  e. 
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Another  slight  modification  has  also  been  made  in  PIPgTEM.  Tn 
the  first  card  of  a  problem  deck,  setting  TI.UIT  =  -500.  causes  the 
program  to  search  for  and  read  initial  temperature  data  as  described 
earlier,  while  setting  TINIT  to  any  value  greater  than  -^85.  causes 
the  program  to  assume  a  constant  initial  temnerRture  eoual  to  the 
value  of  TINIT.  The  new  capability  is  that  setting  TIUJT   =  -'100. 
causes  the  program  to  start  with  initial  temperatures  equal  to  the 
final  temperatures  for  the  preceding  problem.  Unless  the  values  of 
N  are  equal  for  the  two  problems,  some  nonsense  is  generated.  In 
use,  of  course,  this  capability  is  intended  for  a  succession  of  two 
or  more  problems  involving  the  same  pipe,  divided  into  the  same  num- 
ber of  layers.  A  modified  set  of  OTIMENT  cards  at  the  beginning  of 
the  program  describes  the  modified  capabilities  and  the  details  of 
the  problem  deck. 
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